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Abstract 

We determine the abundance of the lightest (dark matter) sterile neutrinos created in the 
Early Universe due to active-sterile neutrino transitions from the thermal plasma. Our start- 
ing point is the field-theoretic formula for the sterile neutrino production rate, derived in our 
previous work [JHEP 06(2006)053], which allows to systematically incorporate all relevant 
effects, and also to analyse various hadronic uncertainties. Our numerical results differ mod- 
erately from previous computations in the literature, and lead to an absolute upper bound on 
the mixing angles of the dark matter sterile neutrino. Comparing this bound with existing 
astrophysical X-ray constraints, we find that the Dodelson-Widrow scenario, which proposes 
sterile neutrinos generated by active-sterile neutrino transitions to be the sole source of dark 
matter, is only possible for sterile neutrino masses lighter than 3.5 keV (6 keV if all hadronic 
uncertainties are pushed in one direction and the most stringent X-ray bounds are relaxed by 
a factor of two). This upper bound may conflict with a lower bound from structure forma- 
tion, but a definitive conclusion necessitates numerical simulations with the non-equilibrium 
momentum distribution function that we derive. If other production mechanisms are also 
operative, no upper bound on the sterile neutrino mass can be established. 
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1. Introduction 

In recent works [TJ H] it has been demonstrated that an extension of the Minimal Stan- 
dard Model (MSM) by three right-handed gauge-singlet fermions — sterile neutrinos — with 
masses smaller than the electroweak scale, allows to address a number of phenomena which 
cannot be explained in the framework of the MSM. Indeed there exists a parameter choice 
within this model, called the i/MSM in refs. [U |2], which is consistent with experimentally 
observed neutrino masses and mixings, provides a candidate for dark matter particles, and 
can explain the baryon asymmetry of the Universe [2] , through a CP-violating redistribution 
of the lepton number among active and sterile flavours [3j , followed by an anomalous conver- 
sion of the lepton number in active flavours into a baryon number [1]. Adding to this model 
one neutral scalar field allows to accommodate inflation [5]. Moreover various astrophysical 
problems may find their explanations [6]. 

In the z^MSM, the role of dark matter is played by the lightest sterile neutrino, with a 
mass in the keV range. This dark matter candidate was proposed a while ago in ref. [7] from 
different considerations, and studied later on in a number of works [H]-[T5]. It was pointed 
out in ref. [7] that the couplings of the dark matter sterile neutrino to charged leptons and 
active neutrinos can be so weak that it never equilibrates in the Early Universe. As has 
been stressed recently |16j-|21j. this means that even if the exact values of the mixing angles 
between sterile and active neutrinos were known, the primordial abundance of the dark matter 
sterile neutrinos cannot be predicted. The reason is trivial: the kinetic equation describing 
sterile neutrino production is a first order differential equation in time, and in order to solve 
it one must specify the initial condition, on which the solution depends significantly, because 
inverse processes are much slower than the rate of Universe expansion. In more physical terms, 
one has to know the physics beyond the z^MSM to set the initial conditions; an example of a 
complete framework, involving the inflaton, can be found in ref. [5]. 

Though a prediction of the primordial abundance of the dark matter sterile neutrinos 
requires physics beyond the i/MSM, an upper limit on the mixing angle 9 between sterile and 
active neutrinos as a function of the sterile neutrino mass M s , 

sm 2 (29) < f(M s ) , (1.1) 

can be established. Indeed, some number of sterile neutrinos are certainly produced in the 
Early Universe through the mixing with active neutrinos, and this amount must be smaller 
than the abundance of dark matter, known from observations. This upper bound can be 
derived with a small number of extra assumptions, which we formulate as follows: 

(i) The i/MSM is a good effective theory below energies of a few GeV, so that there are no 
interactions beyond those included in the z/MSM Lagrangian. 

(ii) The standard Big Bang scenario is valid starting from temperatures above a few GeV. 
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(iii) The charge asymmetries of the plasma (particularly the asymmetries in the total lepton 
number, and in the various lepton flavours) are small, i.e. at most within a few orders 
of magnitude of the observed baryon asymmetry, at temperatures below a few GeV. 

(iv) The masses of the two heavier sterile neutrinos are large enough so that they decay 
above temperatures of a few GeV. 

The assumptions (i— iii) are crucial. For example, if the inflaton is light, its interactions 
may result in sterile neutrino production at a low scale [5], making (i) invalid. For very low 
reheating temperatures after inflation, the assumption (ii) is not satisfied and the production 
of sterile neutrinos may be suppressed [22]. As for point (iii), it is known from ref. [8] that 
the rate of sterile neutrino production is in fact greatly boosted if relatively large lepton 
asymmetries (corresponding to chemical potentials fi/T > 10~ 5 ) are present. On the other 
hand, if the crucial assumptions (i— iii) are valid, the assumption (iv) can be relaxed: once 
the result is known, the decays of heavier sterile neutrinos can be taken into account a 
posteriori [19]. We will return to the last point below. 

There are quite a number of computations of the function f{M s ) already existing in the 
literature, starting from similar assumptions [TJ [TUl ttSJ [13]. They are based, however, 
on kinetic equations which were not derived rigorously and which in fact differ from the 
expressions obtained from the first principles of statistical mechanics and quantum field 
theory [21]. Moreover, as discussed in ref. [21] . neither the hadronic scattering contributions 
to the rate of sterile neutrino production, nor the uncertainties in the hadronic equation-of- 
state, have been exhaustively analyzed in refs. [TOj, [12], [13]. 

The aim of the present paper is therefore to apply the general formalism of ref. [21] to find 
the amount and spectrum of the sterile neutrinos created in active-sterile transitions. Our 
results differ from those that have previously appeared in the literature, even though it turns 
out that the order of magnitude remains the same. 

Given the results of the theoretical computation, avoiding the overclosure of the Uni- 
verse allows us to establish the constraint in Eq. (|l.ip . We can, however, also carry out 
a comparison with other astrophysical and cosmological limits on the properties of sterile 
neutrinos. Indeed, the sterile neutrino radiative decays N — > ^7 produce a feature in the 
diffuse X-ray background [HI [16] or a line in the X-ray spectrum in the direction where dark 
matter is accumulated (such as clusters of galaxies |11[ [T8] , dwarf galaxies [20] , or galaxies 
[11[ [2"0"1 [23| 124]). The position of this line, if found, determines the mass M s of the sterile 
neutrino, while the line intensity would fix the mixing angle 8. No feature or line has been 
observed so far, which places a constraint on the mixing angle of the form sin 2 (26) < fx(M s ) 
[l6l[ll[20l[23l[H[23[Ml[2a[Ml[29]. An exclusion plot from refs. [2Ql [26] is reproduced in 
Fig. [9] below, and can be compared with Eq. (jl.ll) . 

Another, completely independent constraint comes from cosmological structure formation, 
particularly in the form of Lyman-a forest observations [30]— [33] . Being relatively light, 
dark matter sterile neutrinos would play the role of warm dark matter, with a free-streaming 
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length exceeding greatly that of cold dark matter. This erases inhomogeneities on the smallest 
scales. An observation of the small scale structures puts, therefore, an upper bound on the 
free-streaming length and, consequently, on the average velocity of the dark matter particles. 
This converts to a lower bound on the inverse velocity, which can be expressed roughly as 
M s (|q a |)/(|q s |) > Mo, where (|q a |) and (|q s |) are the average momenta of active and sterile 
neutrinos, respectively, at the moment of structure formation, and the value of Mq is as large 
as Mq ~ 14.4 keV according to ref. [32]. The computations of refs. (TJUniE] lead to a nearly- 
thermal spectrum of sterile neutrinos which is somewhat shifted in the infrared in comparison 
with the Fermi-Dirac distribution. According to ref. [13], (|q s |)/(|qa|) — 0.9 (uncertainties in 
this number and its dependence on the sterile neutrino mass will be discussed below). This 
leads to M s > 13 keV [32] . to again be compared with Eq. (jl.ip . Note that ref. [33J gives 
a somewhat weaker bound M s > 10 keV. We stress that these mass bounds do depend on 
the mechanism of sterile neutrino production through the momentum distribution function 
of sterile neutrinos [19]. A model-independent constraint (the Tremaine-Gunn bound) on the 
mass coming from the analysis of the rotational curves of dwarf satellite galaxies |34[ [35] is 
in fact much weaker and reads M s > 0.3 keV [36] . 

These three different types of constraints, Eq. (jl.ip . X-ray, and Lyman-a or the Tremaine- 
Gunn bound, determine the allowed values for the mass and mixing angles of the dark matter 
sterile neutrino, necessary for planning for their search in space missions [29] and in laboratory 
experiments [37]. In addition, the validity of the Dodelson-Widrow scenario, where only 
thermal production is taken into account (so that Eq. (jl.ip becomes an equality), can be 
tested. 

The plan of the paper is the following. We start by reviewing the formalism of ref. [21] 
in Sec. [21 In Sec. [3] we apply it to the leptonic contribution to the sterile neutrino production 
rate, while in Sec. [4] the kinetic equation for the sterile neutrino abundance and its solution 
in the cosmological context are discussed. In Sec. [5] we analyse the hadronic contributions 
that play a role in the solution of the kinetic equation, and their uncertainties. In Sec. [6] 
we combine all these results together and present numerical estimates for the relic sterile 
neutrino abundance as a function of M s and 9. In Sec. [7J we compare our results with the 
observational constraints mentioned above, and discuss the viability of the Dodelson-Widrow 
scenario for sterile neutrino production. We conclude in Sec. In this paper the notations 
of ref. [21] will be used unless stated otherwise. 

2. Review of the general formalism 

As has been discussed in ref. [21], the assumptions (i-iv) allow to derive, from first principles, 
a general formula for the sterile neutrino production rate, provided that we choose parameter 
values consistent with the observational constraints mentioned above. The restriction to 
temperatures below a few GeV in these assumptions means that the electroweak symmetry 
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is broken, whereby it is an excellent approximation to replace the Higgs field by its vacuum 
expectation value v ~ 246 GeV. Moreover combining the two sets of observational constraints 
restricts the mixing parameter(s) 9 to be small. To be precise, there are several mixing 
parameters, 9 a i, where I is the sterile neutrino flavour and a is the active neutrino flavour. We 
define 9 2 y j = \MD\ 2 a j/M] , where |M£>| Q / = \vF a j\/ '\/2; F a i are the neutrino Yukawa couplings; 
and Mi are the sterile neutrino Major ana masses. In the following we choose the lightest 
sterile neutrino to correspond to I = 1. Because of the smallness of 9 a j, 9 a i < 10~ 3 (or, in 
terms of Yukawa couplings, \h a \\ < 10 , \h a 2\, |/i«3 1 ^ 10 -7 ), it is perfectly sufficient to 
restrict to leading order in a Taylor series in 9^j, which simplification plays an essential role 
in the first-principles derivation presented in ref. [21] . (For I = 2,3 these constraints derive 
from the baryon asymmetry of the Universe [2].) For future reference, let us also define 

a=e,fi,r 

this angle corresponds to that in Eq. (|l.ip . 

With these premises, the phase space density ni(t, q) of sterile neutrinos in either spin 
state, 

, . v-^ divj s) (t,x,q) , , 

obeys the equation 

d „ d 



- - Hqi-^-J m(t, q) = R(T, q) , (2.3) 

where H is the Hubble parameter, H = dlna(t)/d£, and qi are the spatial components of q. 
Note that in thermal equilibrium, ni(t, q) = 2nF(g°)/(27r) 3 , where uf is the Fermi distribution 
function. The source term reads 1211 



P (T n) _ 4n F (g°) f \Md&i f . 

1 ,W (2vr) 3 2g°^ 1 {[Q + ReS] 2 -[ImS] 2 } 2 + 4{[Q + ReS]-ImS} 2 { ' 



x 



Tr{^a L (2[Q + ReS] -ImS [Q +Re^] - {[Q + ReS] 2 - [ImS] 2 } lm^)a R ) , 



where Q is the on-shell four-momentum of the sterile neutrino (i.e. Q 2 = Mf); S = T, aa (Q) 
is the self-energy of the active neutrino of flavour a; and aL,o,R are the chiral projectors. 
It is obvious that n/(t,q) and R(T,q) are functions of q = |q| only, and we will use the 
corresponding simplified notation in the following. 

Now, the real part ReS of the active neutrino self-energy is generated at 1-loop level 
through W and Z-boson exchange. At low energies, the result can be expanded in l/m 2 ^, 
l/m?7. In the absence of leptonic chemical potentials, the first term in the expansion vanishes, 
and the leading contribution comes from the second term. Writing the result as 

Re £ ca(Q) = Q a aa (Q) + ib aa {Q) , (2.5) 
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Figure 1: The function b aa (Q) that determines the real part ReS aQ (Q) (cf. Eq. (I2.5|) l. in units of 
G^T 4 q°, as a function of the temperature T and the active neutrino flavour a, with a = e, n, t. We 
have assumed here that |g°| <C raw 



where u = (1, 0), we note that the function a aa {Q) can be ignored, since it is small compared 
with the tree-level term Q . On the other hand the latter structure in Eq. (|2.5h does not 
appear at tree-level, and needs to be kept. For q <C m\v it reads [38l 139] 



b aa (Q) 



16G 2 



7TQ,, 



V 



2(j)(mi a ) + cos 2 6* w (f>{m Va 



(2.6) 



where Gf = g^j^^Flm^ is the Fermi constant, mi a is the mass of the charged lepton of 
generation a (li = e,l2 = = t), and m Ua = is the mass of the MSM active neutrino. 
The function <p is finite and easily evaluated numerically: 



4>{m) 



d 3 p n F (E) 
(2vr) 3 2E 



I 1 2 i 2 

|p| + m 



(2.7) 



Note, in particular, that 0(0) = 77r 2 T 4 /360. The functions b aa are plotted in Fig. [TJ 

As concerns the imaginary part ImS, there is a 1-loop contribution from the same graphs 
as for ReS, but it is exponentially suppressed, ~ exp(— mw/T) (cf. Sec. 13. lh . Therefore, at 
low temperatures, the dominant contribution is generated by those 2-loop graphs which are 
not exponentially suppressed. The dependence on Q is more complicated than in Eq. (12. 6ft : 
in general ImS(Q) = G F T 5 f(Q,T), where / is a non-trivial dimensionless function, which 
is numerically of order unity. 
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Given that Im E <C ReS because of the ct^-suppression, and that Re £ » u b aa (Q) because 
of the reasons mentioned above, the expression in Eq. (|2.4p can be simplified. Carrying out 
the Dirac traces; combining equivalent terms; and re-introducing Dirac traces and chiral 
projectors aL,a>R around ImX at the end, as a reminder of our convention of showing them 
explicitly [2Tj, we arrive at 



(Mf - a L Tm$a R } + 2{Mj + q°b)bTi {fia L Im$a R } 

X } T2 ' ( 2 - 8 ) 

[Mf + 2q% + b 2 ] 

where b = b aa (Q), and we made use of Q 2 = Mf. 

Now, let us estimate the orders of magnitude that are relevant for Eq. (|2.8p (for analogous 
earlier discussions see, e.g., refs. [3(9]). We consider momenta of order |q| ~ T, where the rate 
turns out to peak; temperatures in the range 1 MeV < T < 10 GeV; and sterile neutrino 
masses in the range 10 _1 keV < M 1 < 10 3 keV. Thereby q° ~ T » M u and b ~ 50G 2 F T 5 
(cf. Fig. Q]). We now identify two different temperature scales: 

• We define T x such that M\ ~ 2q% for T ~ T 1} i.e. 

Note that for T ~ Ti, b 2 < M|, and obviously also 6 < g°. 

• We define T 2 such that M\ ~ 6 2 for T ~ T 2 , i.e. 



i 



Note that for T ~ T 2 , » Mf, while it is still true that b < g° 



Let us then estimate the magnitude of the second row in Eq. (|2.8p for various temperatures. 
For simplicity we set ImS ~ G 2 F T^ , since the precise magnitude plays no role, given that 
ImS appears linearly in all terms in the numerator. We observe that: 

• For T < Ti, M 2 dominates in magnitude over all thermally generated terms. The 
second row of Eq. (|2.8p then evaluates to ~ G'pT® /M 2 , i.e. decreases fast at low T. 

• For T ~ Ti, the second row evaluates to ~ 1/10 2 . This is the peak value. 

• For T S> T\, the denominator is dominated by (2q°b) 2 and thus starts to increase. The 
second row of Eq. (US} then evaluates to ~ M?/W 4 G 2 F T 6 , i.e. decreases again fast. 
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• Once T ~ T2, the two terms in the numerator of Eq. f|2.8[) are of the same order of 
magnitude. In fact, the first term becomes negative, but this effect is compensated for 
by the second term, so that the expression as a whole remains positive. The overall 
magnitude is now ~ G 2 F T$, i.e. about 1/10 6 of that at T ~ T\. 

• For T S> T 2 , the second row increases as ~ G F T 4 . Given that Gp « 1/(290 GeV) 2 , 
however, these contributions remain very small in the region T < 10 GeV. Moreover 
in practice our rough approximation turns out to be an overestimate, when compared 
with the exact numerical evaluation. 

To confirm this qualitative picture, we have verified numerically that if one starts the 
integration (to be specified in Sec. 0|) from a high temperature T ~ 10 GeV 3> T2, then 
both terms in the numerator of Eq. f)2.8f) need to be taken into account, in order to obtain 
a positive production rate (in the range T > (5. ..10) GeV the exponentially suppressed 1- 
loop contribution to ImS needs to be taken account as well, but this does not change the 
conclusions). Their combined contribution from the high-temperature region is, however, 
completely negligible. For the dark matter sterile neutrino it is therefore in practice enough 
to start the integration from T ~ T2 ~ a few GeV, and keep only the first term in the 
numerator, which is what we will do in the following. 

3. Leptonic contribution to the production rate 
3.1. 1-loop effects 

We now proceed to determine ImS. We assume that the temperature is well below the 
temperature T ew where the electroweak crossover takes place, T ew ~ 100 — 200 GeV. In 
this situation the Higgs phenomenon provides for a good tool for carrying out perturbative 
computations. The simplest contributions to ImS aQ come from the same 1-loop graphs as 
were considered for ReS aQ , involving a single or Z°-propagator; however we do not 
carry out any expansion in the inverse or Z° masses like for Re Tj aa . A straightforward 
computation in Feynman gauge yields 



Im# 



1-loop 
lep 



(Q) = ira w n 



c=w,z 



I 





Q)riFinB2+ 2~3>-«-q 
Q) n B 2(l - ^fi) + 2 — ^-Q 




(3.1) 
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where n F \ = n F (Ei), ubi = nft{Ei); and pw = %i Pz = cos 2 9 W are the "weights" of the 
charged and neutral current channels. Furthermore, Pj = (Ei,pi) are on-shell four-momenta, 

Ei = Jp\ + rf~ c ' ^ 2 = ^P2 + m c » ( 3 - 2 ) 

where m/^, = m/ a and = m Ua = 0. The graphs in Eq. (|3.3p illustrate the various 
processes, with the wavy line indicating the weak gauge boson. The phase space integrals 
remaining can be evaluated numerically as explained in Appendix A. 

For the masses that we are interested in, mi a ,Mj <C mw-i only one of the channels in 
Eq. (|3.ip gives a non-zero contribution, namely the second one. 

Now, as discussed in the previous section, the dominant production of Nj takes place 
at temperatures much below mw, say, T < my^/10. In this case the phase space factor 
ub2 guarantees that the 1-loop contribution is suppressed by at least ~ exp(— myy/T), and 
thus vanishingly small. One might worry that thermal effects on the boson mass make 
mw(T) smaller than naively expected; at the very small temperatures we are infested in, 
however, this effect is negligible. Therefore, the 1-loop contribution is indeed inessential for 
our purposes, and does not play a role in the following. 



3.2. 2-loop effects 



We now move to specifying the 2-loop contribution to ImS in Eq. (|2.8p . originating from 
intermediate states containing leptons only. These effects can be reliably treated within per- 
turbation theory, and are not exponentially suppressed at low temperatures. At the same 
time, previous evaluations in the literature have made used of phenomenological approxi- 
mations which are not part of the strict perturbative computation, particularly in order to 
simplify the Dirac structures that enter the sterile neutrino production rate. In the following 
we evaluate the leptonic contributions without any such approximations. 

In order to proceed, it is actually helpful to start by recalling the contribution that emerges 
from a pair of free quarks, with masses m 2 , m 3 ; the derivation of this result has been discussed 
in explicit detail in ref. [21]. The final result, given by Eq. (3.41) of ref. [21], reads: 



Im^ad OP (Q) 



(2vr)V 4 
+ (2vr)V 4 
+ (2vr)V 4 
+ (2vr)V 4 
+ (2tt)V 4 



d^pi 

c=w,z 

Q) nFin F2 nF3 A(-mi c , m 2 
Q)n F2 n F3 (l - n F1 )A(mi c ,m 2 



2N c G 2 F n F \ q «) V pc [ - d ' Pl f d " P2 / 
F n=£r J (2vr) 3 2£i J (2tt) 3 2£ 2 J 



d 3 P3 



(2nf2E 2 J (2vr) 3 2£ 3 
m 3 ) + 2-^* q 



Pi + P2 + P3 

P 2 +P3-P1-Q) n F2 n F3 (l - n F1 )A(mi cl m 2 , -m 3 ) + 
Pi + P 3 - P 2 - Q) n Fl n F3 (l - n F2 ) A(-mi c , -m 2 , -m 3 ) + 



Pi +P 2 - P 3 -Q) n F1 n F2 (l - n F3 )A(-mi c ,m 2 ,m 3 ) + 



Pi - P 2 - P 3 -Q) n F i(l - np2)(l - n F3 )A(-mi c ,-m 2 ,m 3 ) + i-*^q 

• 3 
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+ (27r) i 5^(P 2 -P 1 -P 3 -Q)n F2 (l-n F1 )(l-n F3 )A(m lc ,m 2 ,m 3 ) + 2^C-q 
+ {2n) 4 5^\P 3 -P l -P2-Q)n F3 (l-n F1 ){l-n F2 )A(mi c ,-m 2 ,-m 3 )+ s^Q 

+ (2vr) 4 5( 4 )(-P 1 - P 2 - P 3 - Q) (1 - n F i)(l - «F 2 )(1 - n F3 ) «4(m; , -m 2 , m 3 ) }, ^ 

(3*3) 



where npi = n F (Ei); 



A(mi c ,m 2 ,m 3 ) = 7 M (fi + mi c )j u Tr (f 2 + m 2 )7 / ,r(f > 3 + m 3 )7^r 



(3.4) 



Ply = 2, pz = 1/2 are the "weights" of the charged and neutral current channels; and 
mi c has the same meaning as in Eq. (|3.2|) . Furthermore, Pj = (E^p^) are again on-shell 
four-momenta, now with the mass assignments 



Ei = J Pi + m 2 



la 



E 2 = VP2 + m 2 . ^3=\/P3 + m 3- 



(3.5) 



The graphs in Eq. (|3.3f) illustrate the various processes. Further details (in particular the 
values of the masses m 2 ,m 3 and of the Dirac matrix T) will be explained presently. 

Considering then the leptonic contributions, it is obvious that the same eight kinematic 
possibilities will appear as in Eq. (|3.3p . just with different masses and coefficients. On the 
other hand additional terms appear as well, since Z-exchange can proceed through more 
channels than in the hadronic case. 

Given that the kinematic possibilities are identical, however, it is possible to "factorise" 
the result into a "Dirac part" and a "kinematic part". It is therefore enough to show the 
result for the Dirac part by considering one of the kinematic channels only; for simplicity we 
choose the alternative in Eq. (|3.3|) where all masses are positive inside the ^4-function. A 
straightforward computation parallelling the one in ref. [21] then produces the result 



Im?j^(Q)=2GSV( 9 ' 



d 3 P i 



d 3 P2 



d 3 P3 



(2vr) 3 2£i J (2tt) 3 2E 2 J (2tt) 3 2E 3 



x | 2 7 ^i + ro, a )7" E ^[{f2 + mi fj )l»a L {P 3 + m V p) lv a L 

P=e,n,r 

+ l^(P 1+r n v JY E TiUp 2 + m h )^T(P 3 + m h h u r 



+ 



+ 



(3=e,fi,r 



+ ^7 M (f l + m Va )Y E Tr [(^2 + m U0 ) ltl a L (P 3 + m u ^ u a L 

/3=e,n,T 



- -f{f i + m Va ) 1 u a L {P 2 + m Va ) 1 ^a L {P 3 + m Va ) lv 
r(Pi + m la )Y r (f 2 + m la )^a L (P 3 + m Va ) lv - 
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x(27r)V 4 )(P 2 -Pi-P3-Q)n F2 (l-n F i)(l-n F3 ) + ... , a-^Q (3.6) 

^ 3 

where now F = —1/2 + 2x w + 75/2, with x w = sin 2 W , and the energies E{ are on-shell with 
the obvious mass assignments. 

Inspecting Eqs. f|2.8j) and (|3.6p . we observe that two kinds of Dirac traces appear in the 
final result: either a product of two traces, or a single trace. Both cases are elementary, and 
result in 



Ti = TV 



Tr 



T 2 



16 
Tr 



2 + m 2 )7 M (a + 675)^3 + "13)7^(0 + ^75) 
(a - 6) 2 Pi • P 2 P 3 • E + (a + bf P 1 ■ P 3 P 2 • P + (6 2 - a 2 ) m 2 m 3 Pi • p] , (3.7) 

Ea Ll »{P l + mi)7"(o + 67s) 2 + m 2 )^(c + d 75 )(P 3 + 7713)7^ 
-2(a + 6)(c + d) Pi • P 3 P 2 • P + (a + 6)(c - d) m 2 m 3 Pi • E + 



+(a — b)(c — d) mim 3 P2 • P + (a — 6)(c + d) m\m2 P 3 • E 



(3.8) 



where the "external" four- vector P is either Q or u (cf. Eq. (|2.8p ). 

Given these ingredients — the kinematic channels in Eq. (|3.3p . the Dirac structures in 
Eq. (|3.6p . and the traces in Eqs. (|3.7p . (|3.8p — we can finally collect together the full expres- 
sion needed in Eq. (|2.8|) . We obtain 



Tr 



d 3 Pi 



d 3 P2 



(2vr) 3 2Pi J (2vr) 3 2P 2 7 (2vr) 3 2P : 



d 3 P3 





(2vr)V 4 ) 


(Pi + P 2 + P 3 - 


+ 


(2vr)V 4 ) 


(P2 +P3- 


Pi ~ 


+ 


(2ir) 4 6& 


(Pi +P3- 


P2- 


+ 


(2vr)V 4 ) 


(Pi +P2- 


P 3 - 


+ 


(2n) 4 S^ 


(Pi - P 2 - 


P3- 


+ 


(2vr)V 4 ) 


(P 2 ~ Pi ~ 


P3- 


+ 


(2vr)V 4 ) 


(P 3 ~ Pi ~ 


P2- 


+ 


(2vr)V 4 ) 


(-Pi ~ P2 


~ Pc 



3 U "- q 



■ 2 

Q 

• 3 



Q 

• 3 



(3.9) 
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where % equals T\ or T2, as specified in Table 1. The prefactors Cj and the masses mi, m2, 
7713 relevant for each channel are also listed in Table 1. For completeness and future reference, 
we have included in this Table the perturbative hadronic contributions as well. 



4. Solution of the kinetic equation 



Now that the ingredients entering Eq. (|2.3p are in place, we need to discuss the solution of this 
equation, and identify the precise information about the equation-of-state which is needed 
for the solution. The considerations that follow are rather standard (see, e.g., refs. [7J [9]), 
but we present them here for completeness and in order to fix the notation. 

It is convenient to represent the rate R(T, q) as a function F of dimensionless variables, 

R(T, q ) = G 2 F T 5 F(^l) , (4.1) 



T T , 

where rrii are the masses of the particles of the MSM and of the sterile neutrinos. We also 
introduce the effective numbers of massless bosonic degrees of freedom g e s{T) and h e d(T) via 
the relations 

tt 2 T 4 2tt 2 T 3 
<T) = -go-WfCT) , S (T)e— MT), (4.2) 

where e(T) and s(T) are the energy and entropy densities, respectively. Given the equation- 
of-state of the plasma [i.e. the relation between the pressure and the temperature, p = p(T)], 
9eff(T) and h c g(T) can be found from the standard thermodynamical relations 

,m\ 30 d fP\ , 45 dp , A n . 

9^(T) = ' MT) = 2^T3 dT • (43) 

Furthermore the sound speed squared, c 2 s (T) = p'(T)/e'(T) = p'(T)/Tp"(T), can be expressed 
as 

d{T) 6+ h cS (T) ' (4A) 
where a prime denotes a derivative with respect to T. The Hubble rate is given by 



2 - 1 



H = -f— , M (T) = M PI 



45 



(4.5) 



M (T) ' > ^ l^ 9cS (T) 

where Mp\ is the Planck mass. Note that, from the assumption (iv) in Sec. [H the heavier 
sterile neutrinos do not contribute to g e ff and h e g for the temperatures of interest. Further- 
more, the contributions from the lightest sterile neutrino are negligibly small as long as it 
does not get thermalized. 

The kinetic equation (|2.3p can easily be integrated by writing ni(t,q) = f(t,y), where the 
variable y = a(t)q accounts for red-shift. After simple manipulations one gets the distribution 
function of sterile neutrinos at temperature To: 

G\ Z" 00 , Mq(T)T 2 



3 J t cj(T) 



fm q 


r h eS (T) - 


1 




h e s(T )_ 





(4.6) 
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Table 1: The coefficients and masses that appear in Eq. (|3.9[) . in the perturbative limit. The functions 
7i and Tj, which depend on the masses mi, 771,2, "73 and the coefficients a, b, c, d (which take values 
as specified above), are defined in Eqs. (13. 7| and (|3.8|) . respectively. The symbols / a stand for leptons 
of generation a, i.e. m; 1 = m e , 777; 2 = m M , m; 3 = m T . 
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The integral of Eq. (|4.6p over the momenta q gives the number density of sterile neutri- 
nos, which can conveniently be normalized with respect to its would-be equilibrium value, 
n eq (io) = 3C(3)T 3 /2^ 2 : 

_ 8GX f°° dzz 2 M (T)T 2 hes(To) F /ro^ \ _ 



n eq (t ) 9C(3) y To Jo c|(T) h cS (T) \T 

Finally another characterization of the same quantity is obtained by normalizing through the 
total entropy density, which produces the so-called yield parameter: 

s(to) 7T Jt Jo cj(T)h eS (T) \T J 

The benefit of Eq. (|4.8p is that it obtains, unlike Eq. (14. 7h . a constant value at low tempera- 
tures (in the absence of entropy production), because the factor /i e ff(To) drops out. 

To summarise, to find the relic concentration of sterile neutrinos produced by active-sterile 
transitions, one should have a reasonable approximation for the hadronic equation-of-state 
p(T) and for its first and second temperature derivatives, p'(T) and p"(T). 

To convert the result of Eq. (j4.8f) to physical units, we denote the contribution to the 
number density of the lightest sterile neutrino (1 = 1) from an active neutrino of flavour a, 
by n a i (so that n\ = J2a=e,fj,,T n ai)- To relate this quantity to Qdm, we write 

O al . = ^ , (4.9) 

Per Per/ S 

where Y a \ = n a \j s. From Particle Data Group [51], one finds p cr ~ 1.054 x 10 -5 h 2 GeVcm~ 3 
and s = 7.04n 7 ~ 2886cm -3 , yielding p C r/s ~ 3.65 x 10 -9 /i 2 GeV . Then the result of 
Eq. (14. 8p can be expressed as 



!i«i/. 2 = 0.11C t ,(A ! f 1 )( ! |^f-) , Mid) 



2 



where C a (Mi) = 2.49 x 10~ 5 x Y al /6 2 al x (keV/Mi). 

The total relic density of the lightest sterile neutrino is now given by Y^oc=e,u,r ^aih 2 . On 
the other hand, the dark-matter density in the present Universe has been presicely measured 
by the WMAP collaboration [40J (68% CL): 

n dm h 2 = 0.1051^13 . (4.11) 

Therefore, to avoid the overclosure of the Universe by N\, we get Y^a=e,u,T Qaih? ^ 0.112, 
which leads to 

/|M D | al x " 



This corresponds to Eq. (|l.ip . In the following, we will present our results for the functions 
C a {M x ) defined by Eq. flUU]). 



13 



100 



0.4 




r/MeV T/MeV 



Fi gure 2: Left: <7effj^eff a s defined in Eq. (|4.3p . for the MSM and for the QCD part thereof, for 
T = 1 MeV ... 4 GeV (for more details, see ref. [43]). Right: the speed of sound squared c 2 s , for the 
same systems. Various sources of uncertainties in these estimates are discussed in the text. 

5. Hadronic uncertainties 

As discussed in ref. [21 j . there are in principle three sources of hadronic effects and cor- 
responding hadronic uncertainties in the present computation: those entering ReS, ImS, 
as well as the overall equation-of-state through the considerations in Sec. 01 The hadronic 
contributions to ReS are suppressed by a w with respect to the leptonic ones and will be 
omitted in the following. The hadronic contributions to ImS arise at the same order as 
the leptonic ones, consituting thus a significant source term. Even more important are the 
hadronic contributions to the equation-of-state, since in the temperature range of interest 
hadrons completely dominate over leptons as far as the change in the effective number of 
light degrees of freedom is concerned (cf. Fig. [2|). On the other hand, the hadronic effects in 
the equation-of-state are understood somewhat better than the ones in ImS, given that they 
are theoretically more straightforward to access, both in the framework of the weak-coupling 
expansion [41j-|43j and of lattice simulations [44j-[47j, so that the relative uncertainties are 
perhaps smaller then in ImS. In fact we estimate the absolute uncertainties from these two 
sources to be of similar magnitudes. In the following we discuss the ways in which we extract 
these two quantities and try to estimate their uncertainties. 
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5.1. QCD equation-of-state 

A basic fact to realise about the QCD equation-of-state is that the low-temperature hadronic 
world and the high-temperature partonic world can be analytically connected to each other, 
by tuning quark masses and the temperature: at non-zero quark masses there is no order 
parameter to distinguish between the two situations. Moreover, lattice simulations have been 
suggesting for quite a while already that such an analytical "crossover" is met even if the quark 
masses are not tuned but kept at their physical values; for a recent studies see, e.g., refs. [38] , 
In fact, these statements could even have been made without recourse to lattice simulations: 
the thermodynamics of a dilute "gas" of hadronic resonances at low temperatures, and of 
interacting quarks and gluons at high temperatures, appear to extrapolate to each other 
surprisingly well, when determined precisely enough [43J. 

Now, even if there only were a smooth crossover in the system, it is conventional to assign 
a (pseudo) critical temperature to it; we will denote this temperature by T c . Of course T c 
cannot be defined precisely, but the "conceptual" ambiguities involved do not appear to be 
much larger than the current statistical uncertainties. A recent large-scale lattice study, for 
instance, suggests that T c ~ 192 ± 8 MeV 09]. 

On the other hand, it is also important to realise that at the current moment the lattice 
studies still involve systematic uncertainties (related for instance to finite-volume effects or 
the absence of a precise continuum extrapolation) , which are not quantitatively under control. 
Consequently, even though the final word on the value of T c and on what happens around 
it lies with the lattice simulations, they are still rather far from establishing the correct be- 
haviour in a wide temperature interval. In fact, even the results of various groups employing 
similar techniques differ by much more than the statistical uncertainties cited above [50] . 
Moreover, lattice studies fail to show any signs of the characteristic peak that physical pions 
cause in the sound speed at T 70 MeV [33] (see also Fig. [2]), a problem which can probably 
be assigned to the fact that the quark discretizations used in the current finite-temperature 
simulations do not respect the chiral symmetry that plays an essential role at low temper- 
atures, and also tend to employ unphysically heavy quark masses. Lattice simulations also 
cannot be applied at very high temperatures. To summarise, the current numbers can be 
expected to be fairly reliable in the range from about T c to about (2...3)T C at best. 

For these reasons, we prefer to adhere to the procedure introduced in ref. [43] here, rather 
than to lattice simulations. It makes use of a gas of hadronic resonances at low tempera- 
tures; the most advanced (up to resummed 4-loop level [H]) weak-coupling results at high 
temperatures; and an interpolation thereof at intermediate temperatures^ Remarkably, the 
temperature interval where an interpolating function is needed in order to sew together the 
two asymptotic functions is fairly narrow, not more than 10 — 20 MeV, and centered around 
T ~ 200 MeV. We will refer to this temperature in the following as T c , but stress that, de- 
spite the curiously good agreement with the crossover temperature suggested by the lattice 

1 We have corrected a minor error in the numerical results of ref. [43] , 
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study mentioned [49J, T c does not need to coincide with any specific definition of T c , given 
the inherent ambiguity in the location of a crossover. The significant benefit of this recipe, 
compared with lattice simulations, is that the results can be evaluated also at arbitrarily low 
and high temperatures without problems. 

Now, this recipe is naturally not exact in the intermediate temperature range; the results 
can be improved when future lattice simulations get closer to the infinite- volume, continuum, 
and chiral limits. In order to estimate the uncertainties in the present implementation, we 
have considered two types of variations of the basic setup: 

• The shape and "sharpness" of the crossover can be changed by adding artificially heavy 
hadronic resonances to the gas, or removing existing ones. 

• The overall location of the pseudocritical temperature can be changed most simply by 
just rescaling the temperature units by some percentual amount. 

It is naturally not difficult to come up with other possible variations as well, but as such 
recipes are purely phenomenological in any case, it is sufficient for our purposes to restrict 
to these most straightforward possibilities. 

We find that the first of these variations has a relatively minor impact on the results. On 
the contrary, a change in T c would lead to visible changes. Qualitatively, a low T c implies 

1 /2 

a low abundance, since the results are inversely proportional to g e g h e g (cf. Eq. (|4.8p ). and 

1/2 

the production then peaks on the partonic side where g e g h e g is higher. Consequently, we 
will assume in the following that the uncertainties of the equation-of-state can be sufficiently 
estimated by a rescaling of T c by 20% in both directions, T c = (160. ..240) MeV, which should 
be a conservative estimate. 

5.2. Vector and axial current spectral functions 

As has been discussed in ref. [21] . the hadronic contributions to ImS can be expressed through 
the spectral functions related to vector and axial currents with various flavour structures. 
The current status of the in general very demanding challenge of reliably determining these 
quantities was also reviewed in ref. [21j . 

Irrespective of all details, however, the general pattern we expect for the hadronic effects 
should be clear: at high temperatures T > T c , hadrons should be reasonably represented 
by free quarks (cf. Eq. (|3.3p ). whereas once the temperature is lowered, confinement sets in, 
and all hadronic effects eventually switch off, given that all mesons and baryons are massive. 
For T <C T c , for instance, the hadronic spectral functions can be determined with chiral 
perturbation theory. It is easy to see that to leading order in the chiral expansion and 
assuming a mass-degenerate limit, only the flavour non-singlet axial current, denoted by A® 
in ref. [21], contributes. We have verified explicitly that this dominant pionic contribution 
is strongly suppressed with respect to the leptonic contributions at all relevant temperatures 
(below 1%), for the small masses Mi that we are interested in. 
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To give another argument in the same direction, we may reasonably assume that the 
hadronic effects are proportional to the number of effective hadronic degrees of freedom, 
/i^ D , related to the hadronic contribution to the entropy, which also rapidly decreases with 
the temperature (cf. Fig. [2]). 

These considerations suggest that a strict upper bound for the hadronic contribution to 
ImS can be obtained by simply computing the effect of free quarks, Eq. We take into 

account the u, d, s and c quarks, cf. Table 1, with their MS scheme masses |51| . A strict lower 
bound can obviously be obtained by simply setting iV c = in the hadronic contributions 
listed in Table 1. As we have already mentioned, the "error bars" that result from this 
(clearly most conservative) procedure are roughly of the same order as those related to the 
equation-of-state, as defined in Sec. 15. 11 (In fact they are a bit larger as we will see, but 
the recipe is also overly conservative.) Furthermore we define a phenomenological "mean" 
value by making use of the perturbative contribution in Table 1, but rescaling N c by the 
factor /iB^ D (T) /58, where 58 counts the hadronic degrees of freedom in the non- interacting 
four-flavour limit. 

6. Numerical results for the sterile neutrino abundance 

For the numerical evaluation of the sterile neutrino abundance it is useful to note that the 
computationally most demanding part entering our equations, ImS, is practically indepen- 
dent of Mi in the range Mi <C T that we are interested in. Therefore, Im E can be evaluated 
once and for all by using any fixed M\. In contrast, the "prefactor" (i.e. the parts of Eq. (12, 8ft 
other than ImS) does depend strongly on the precise value of Mi, yielding a non-trivial de- 
pendence on Mi for the function C Q (Mi) of Eq. (|4.10p . 

As another comment, we note that for q° > there are seven channels in Eq. (|3.9p that can 
be realised. It can be verified numerically, however, that in practice only the 2 — > 2 processes 
give a significant contribution; the other cases have a restricted phase space and amount to 
at most 5% of the 2^2 channels (for typical parameter values in fact much less than 5%), 
which is smaller than our uncertainties. To decrease the numerical cost, we have therefore 
omitted these contributions in the following. It may be noted, though, that these additional 
channels systematically increase the sterile neutrino abundance. 

The way we have numerically evaluated the phase space integrations entering Im X at the 
2-loop level is explained in some detail in Appendix B. Once an efficient implementation is 
available, the remaining integrations (over T, q; cf. Eq. (|4.8H ) are relatively straightforward. 
We note that the g-integration converges exponentially fast, and can in practice be limited 
to values q/T < 10. ..15. 

In Fig. [3] we show an example of the To-evolution of the function C a (Mi), as defined by 
Eq. (|4.10p . We observe that most of the sterile neutrino abundance is indeed generated at 
temperatures of a few hundred MeV for Mi in the keV range. 
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Figure 3: An example of the Tb-evolution of C a (cf. Eqs. (|4.8|) - (|4.10|l ). for a fixed a — e and 
Mi = 10 keV. Shown are the two sources of hadronic uncertainties that are discussed in the text: 
from the equation-of-state (EOS) and from hadronic scatterings. 



In Fig. U] the final low-temperature values of C a (Mi) are shown as a function of a and 
Mi. Importantly, we observe that the dependence on Mi is fairly mild for a wide range of 
Ml, say Mi = (0.1-10 3 ) keV, once expressed as in Eq. (|4.10p . This means, roughly speaking, 
that the relic density of N\ is determined by the Dirac neutrino masses |Md|q.i. According 
to Eq. (|4.12p the Dirac neutrino masses should typically be |Md| q i < 0.1 eV. 

Furthermore, due to the differences between C a , the relic abundance of N% depends on the 
flavour structure of the Dirac neutrino masses |Mrj| a i. We find the (moderate) hierarchy 
C e > > C T . Thus, the largest abundance is obtained when 

\M D \ el ^0, \M D \^ = \M D \ T i =0. (6.1) 

In this case (which we call "case 1"), we obtain the most stringent upper bound on 2 from 
Eq. (I4.12D . On the other hand, the case (which we call "case 2") when 

\M D \ Tl ^0, \M D \ el = \M D \^ = , (6.2) 

gives the smallest relic abundance and the weakest upper bound on 9 2 . 

In Fig. we show the upper bounds on the mixing angle for the above two cases. We also 
compare with the most recent previous computation from the literature [13]. It can be seen 
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Figure 4: The functions C a (Mi), together with estimated hadronic uncertainties, as described in the 
text (see also Fig. [3j. Left: a — e,fi, and only the uncertainties for a — e are shown. Right: a = r. 



that the bound obtained is rather insensitive to the flavour structure of the Dirac neutrino 
masses |Mrj| a i, at least when plotted on a logarithmic scale. 

If we ignore the dependence of C a on Mi and a in Fig. HJ setting C a ~ 0.5, we obtain a 
very simple but useful approximate bound, 

sin 2 (2^4 £ 91, < 8xW^|y\ (6.3) 

More precise expressions, based on fitting the numerical data, will be given in Eqs. (|7.ip of the 
next section. In Fig. [5] the units on the y-axis have been so chosen that a direct comparison 
with the approximate formula in Eq. (j6.3|) is possible. 

An alternative normalization of the sterile neutrino abundance is given in Eq. (14. 7p . and 
the corresponding results are plotted in Fig. It is observed that the abundance generated 
is typically much below its equilibrium value, as expected [recall that in order to avoid 
overclosure, (|M D | a i/0.1 eV) 2 < a few, cf. Eq. <KT2h and Fig.g|. 

Finally, in Fig. [7] we show the momentum distribution, n a i(T,q), for a few examples, 
normalised to the massless equilibrium case, n eq (to,q) = 2nF(<?)/(27r) 3 . As noted in the 
literature |13j , the momentum distribution is shifted towards the infrared compared with the 
equilibrium distribution. For heavy M\ the shift is fairly substantial (note the logarithmic 
scale). Structure formation simulations necessitate the momentum distribution function of 
dark matter as input, and curves such as the ones in Fig. [7] could be used for this purpose. 
However, as mentioned in Sec. [H one may get a rough first impression of the results already 
by just rescaling the average momentum to a shifted value. To allow for such an estimate, 
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Figure 5: The parameter values that, according to our theoretical computation (cf. Fig. 0] and 
Eq. (|4.10[) h lead to the correct dark matter abundance in the Dodelson-Widrow scenario; if additional 
sources are present, sin 2 29 must lie below the curves shown. The grey region between case 1 (lower 
solid line) and case 2 (upper solid line) corresponds to different patterns of the active-sterile mixing 
angles, cf. Eqs. (|6.1|) . (|6.2[) . The absolute upper and lower bounds correspond to one of these limiting 
patterns with simultaneously the uncertainties from the EOS and from hadronic scatterings set to 
their maximal values. The yellow band indicates the result in Eq. (7) of ref. [13] , 2 where we have 
inserted fid m = 0.20, and varied the parameter T c in the range T c = (150. ..200) MeV. 

Fig. displays the average momentum of the sterile neutrinos generated, in comparison with 
the equilibrium value for active neutrinos. 

7. The Dodelson-Widrow scenario 

We are now in the position to discuss the Dodelson-Widrow (DW) scenario for sterile neutrino 
dark matter, which is based on the assumptions (i)-(iv) we stated in Sec. [TJ In this scenario, 
the initial abundance of the dark-matter sterile neutrinos was zero at T > 1 GeV, and they 
were only produced in active-sterile neutrino transitions from the thermal plasma. Thus, in 
terms of Eq. (jl.ip . it is the maximal value of the mixing angle which realizes this scenario, 

2 Ref. ;F3] does not state explicitly the sterile neutrino mass range in which its Eq. (7) should be valid. 
After the eprint version of our paper had appeared, we were informed by K. Abazajian that Eq. (7) of ref. [13] 
was meant to be valid for 0.5 keV < Mi < 10 keV. 
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Figure 6: The sterile neutrino abundance normalised to its equilibrium abundance, n a i/n e q, at To = 
1 MeV. Hadronic uncertainties have the same meaning as in Fig.[3j The combination (|Md| q i/0.1 eV) 2 
has been factored out in analogy with Eq. (|4.10[) and, for better visibility, the results have been 
multiplied by Mi/keV. Left: a = e, fi, and only the uncertainties for a = e are shown. Right: a = t. 



sin 2 (20) = /(Mi). We denote M s = Mi in this case. This scenario can be confronted with a 
number of cosmological and astrophysical observations. 

In order to carry out a comparison with the various X-ray constaints that exist, we replot 
the curves of Fig. [5] in Fig. [9] (for a subrange of mass values), and compare with the regions 
excluded according to refs. [20], [26]. To facilitate numerical estimates, we also note that in 
this mass range (and on the resolution of the logarithmic scale), the theoretical curves can 
be roughly approximated by straight lines, 



log lo [sin 2 20] 

log lo [sin 2 20] ~ 

log lo [sin 2 20] ~ 

log 10 [sin 2 26>] ~ 



-6.9267 - 1.7350 log 10 [M s /keV] 
-7.0862 - 1.8451 log 10 [M s /keV] 
-7.2491 - 1.8359 log 10 [M s /keV] 
-7.6195 - 1.7428 log 10 [M s /keV] 



(absolute upper bound) 
(case 2, mean) 
(case 1, mean) 
(absolute lower bound) 



(7.1) 



These fits improve on the first approximation in Eq. (|6.3|) . The lower bound of the region 
excluded by X-ray constaints will be denoted by fx(M s ); for this data a straight line is a 
worse approximation, but restricting to the mass range M s = (2. ..10) keV, the function can 
be roughly approximated by a parabola, 



log 10 [sin 2 2(9] ~ -3.9499 -9.7859 log 10 [M s /keV] + 3.4902 log 2 [M s /keV] . 



(7.2) 



We then note that the intersection of f(M s ) and fx(M s ) gives an upper bound on M s . It 
is found from Eqs. (|7.ip . (|7.2H that the bound M s < 3.5 keV is obtained by using the mean 
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Figure 7: The distribution functions n a i(to,q), for To = 1 MeV, normalised to the massless equilib- 
rium case, n e q(to,q) — 2np(g)/(27r) 3 . The combination (| A/ z? |ai/0.1 cV) 2 has been factored out in 
analogy with Eq. (|4.10[) . Left: a = e. Right: a = t. These results can be compared with Fig. 1 of 
ref. [13], showing the same distribution functions for 0.3 keV < Mi < 140 keV; although similar at 
first sight, there are significant differences on closer inspection, for instance our curves are monotonic 
functions of q/T even for small M 1; unlike those in ref. [15] . 



predictions for the cases 1 and 2 discussed in Eqs. (|6.ip . (16. 2D of the previous section. For the 
most conservative case with respect to the hadronic uncertainties, the bound can be as large 
as M s < 4.3 keV. Adopting the absolute lower bound on f(M s ) with the maximal hadronic 
uncertainties, and also assuming that the X-ray constraints are still uncertain by a factor of 
two, we set M s < 6 keV as the most conservative limit. As we will show below, possible 
entropy dilution through the decays of the heavier sterile neutrino(s) makes this bound more 
severe, and, therefore, this bound is robust. 

The fate of this scenario now depends on constraints from structure formation simulations, 
making use of Lyman-a observations, which give a lower bound on the mass of the dark matter 
sterile neutrino as M s > (|q s |)/(|q a |)M , where M « 14.4 keV (with (|q s |)/(|q a |) ^ 0.9) [32] 
and M 10 keV (with (|q s |)/(|q a |) ^ 1-0) [33j at 95% CL. Note, however, that these bounds 
are obtained by assuming that the sterile neutrino momentum distribution function is propor- 
tional to the Fermi-Dirac one (with a certain rescaling of the average momentum in ref. |32j), 
even though the sterile neutrino is not thermalised. As shown in Fig. [71 the momentum 
distribution function is in fact significantly distorted from the Fermi-Dirac shape towards the 
infrared (apart from having a smaller normalization). Assuming that the corresponding ef- 
fects can be captured by the factor ~ (|qs|)/(|q a |), we read from Fig.[8]that (|q s |)/(|q a |) — 0.8 
for M s ~ 10 keV. Given that refs. [32], [33] assumed (|qs|)/(|q a |) — 0.9,1.0, respectively, the 
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Figure 8: The average sterile neutrino momentum compared with the active neutrino equilibrium 
value, (|qs|)/(|qa|), where (|q a |) = 7tt 4 Tq/ 180£(3) ss 3.15To. Left: a = e,[i, and only the uncertainties 
for a — e are shown. Right: a = t. We have assumed for this figure that only one of the mixing 
angles 9 a \ is non- vanishing. 



results for lower bounds can be re-interpreted as M s > 11.6 keV and M s > 8 keV. Thus the 
X-ray and Lyman-a regions do not overlap, and astrophysical and cosmological constraints 
appear to rule out the DW scenario, in spite of all theoretical uncertainties involved in the 
computation. The scenario can only survive if the computations of refs. |32[ [33] turn out to 
be uncertain for some reason, be it on the observational side, or on the theoretical side, re- 
lated to the momentum distribution of the sterile neutrino (i.e., the need to use of the proper 
distribution function from Fig. [7] rather than the Fermi-Dirac one). We should mention that, 
even in such a situation, the Tremaine-Gunn lower bound M s > 0.3 keV for fermionic dark 
matter does remain valid 1361 . 



7.1. Entropy dilution 

In the i/MSM, there are two heavier sterile neutrinos in addition to the lightest dark-matter 
one. As shown in ref. [TS], the heavier sterile neutrino(s) could momentarily dominate the 
energy density of the Universe, and their subsequent decay(s) could cause a significant entropy 
dilution after the dark matter sterile neutrinos have been produced. Here we briefly reiterate 
the corresponding effects, i.e., relax the assumption (iv) of Sec. [TJ 

Due to the entropy release, the yield Y\ of the dark matter sterile neutrino is diluted by 
a factor S compared with that without the entropy dilution. Thus, the upper bound on the 
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Mj / keV 



Figure 9: The central region of Fig. [5] (Mi/keV = 0.3. ...50.0), compared with regions excluded by 
various X-ray constraints. From left to right, the X-ray constraints come from refs. [2D], [2DJ, [2D] , 
respectively. The abbreviation LMC stands for the Large Magellanic Cloud, and the names in the 
parentheses refer to the X-ray observatories whose data were employed for deriving the constraints. 

active-sterile neutrino mixing angle becomes weaker, 

sin 2 (2#) < f(M s )xS. (7.3) 

Note that fx{M s ) from the X-ray constraints does not change, since fx is derived from the 
present dark matter density f^dm- Then, we can see from Fig. [9] that the upper bound on M s 
becomes more stringent for a larger S. 

On the other hand, the lower bound on M s from the Lyman-a observations also decreases 
if S > 1, by a factor S" 1//3 if N\ is out-of-equilibrium [19]. The combination of the two effects 
might open a window for sterile neutrino dark matter even in the presence of the Lyman- 
a constraints. To quantify the effect, we note that in the region M s < 4 keV, the X-ray 
bound can be approximated as log 10 [sin 2 26>] < —4.6121 — 6.4659 log 10 [M s /keV]; the most 
conservative bound from Eq. (|7.ip becomes log lo [sm 2 20] < —7.6195—1.7428 log 10 [M s /keV] + 
log 10 [5']; and the Lyman-a bounds obtain the form M s x S 1 ^ 3 > (8. ..11. 6) keV. 

If M s x S 1 / 3 > 8 keV from ref. [33], we find that there appears an allowed region for 
S > 155. When S ss 155, a sterile neutrino with M s ps 1.5 keV and sin 2 (26») ps 1.9 x 10" 6 
could function as dark matter and be consistent with X-ray and Lyman-a constraints. On the 
other hand, if M s x S 1/3 > 11.6 keV from ref. [32], there is an allowed region for S > 3.3x10 , 
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corresponding to M s < 0.78 keV and sin 2 (2#) > 1.2 x 10~ 4 . However, the combination 
(|Mo| a i/0.1 eV) 2 = 2.5xl0 7 (M s /keV) 2 sin 2 (26>) is larger than 10 3 in this case and, according 
to Fig. EJ N\ are thermalized, whereby the assumptions we have made are no longer self- 
consistent. We discuss the correct procedure for the latter case in Sec. 17.21 

To summarise, we have learned that the Dodelson-Widrow scenario remains a possibility if 
there exists a large enough entropy dilution and if the lower bound on M s from the Lyman- 
a data is small enough. However, it is very difficult to get such a large dilution factor 
within the z^MSM. It has been shown in ref. [T9] that the z/MSM induces S ~ 30 at most 
and S > 100 is obtained only if there exists some physics beyond the z^MSM. Therefore, to 
realize this scenario, the constraints from the Lyman-a data should again be weaker than 
those presented in refs. [32j 133]. 

7.2. The case of thermalised sterile neutrinos 

So far, we have assumed that the Dirac neutrino masses |Afo|al are sufficiently small so that 
N\ was never in thermal equilibrium. Now, let us discuss the possibility of thermalized N\ 
as dark matter. If JVi is thermalized at some high temperatures and decouples at T = T^i, 
the present relic density is given by 



Note that this result is independent of 9 as long as 9 is sufficiently large such that N\ was 
thermalized at T > T^i- To explain all of dark matter by such a density of N\s, the sterile 
neutrino mass should be 



The ratio /i e ff/10.75 is at least unity, cf. Fig. [2j In any case, this value is much smaller than 
the Tremaine-Gunn bound, M s > 0.3 keV, and is also below the bounds from the Lyman-a 
constraints for the thermalised case; M s > 2.4 keV [32] and M s > 2 keV [33] (at 95% CL). 

If there exists entropy dilution, the mass of the once-thermalized dark matter sterile neu- 
trino becomes 



On the other hand, the Tremaine-Gunn bound and the Lyman-a constraint remain un- 
changed. This is because the distribution function of the sterile neutrino was the Fermi-Dirac 
one at T = T^i and it is just red-shifted to the present time. Although the Tremaine-Gunn 
bound can be satisfied with S ~ 30, larger factors S > 240 and 200 are required to avoid 
the Lyman-a bounds from refs. [32\ I33j. respectively. Therefore, the thermalized N\ can be 
dark matter within the z^MSM only if there is the largest possible entropy dilution and if the 
lower bound on M s is given by the Tremaine-Gunn bound, rather than the Lyman-a bounds. 




(7.4) 




(7.5) 




(7.6) 
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7.3. Summary 

We have shown in this section that if the lower bound on M s from the Lyman-a data is 
indeed as found in refs. [32^ [33], then there is no parameter space left within the z^MSM to 
realize the scenario of sterile neutrinos serving as warm dark matter, respecting the "crucial" 
assumptions (i)-(iii) of Sec. [TJ 

Of course, this does not mean that the sterile neutrino is excluded as a dark matter candi- 
date, but that at least one of the assumptions would have to be relaxed for this to be the case. 
For instance, sterile neutrinos could be produced mainly in interactions beyond the z^MSM, as 
in ref. [5]. Another logical possibility is that the Universe had substantial leptonic asymme- 
tries at small temperatures, leading to a resonant production of sterile neutrinos [8]. In any 
event, for M s > 3.5 keV (6 keV if all hadronic uncertainties are pushed in one direction and 
the strongest X-ray bounds are relaxed by a factor of two) the active-sterile transitions from 
the thermal charge-symmetric plasma cannot produce cosmologically interesting amounts of 
sterile neutrino dark matter. 

Finally, we should like to stress that, within the corner of the parameter space that defines 
the z/MSM, the bound on the mixing angle sin 2 (2#) < f(M s ), i.e. Eq. (|4.12|) . gives a prediction 
on the mass scale of the lightest active neutrino, m Ul , whose see-saw formula includes the 
lightest Majorana mass M\ only pQ. By setting C a {M\) ~ 0.5, we find 



If N\ is the dark matter particle, its mass should in any case be M\ > 0.3 keV from the 
Tremaine-Gunn bound. Then, we get the weakest bound as m ui < 6.7 x 10~ 5 eV, which is 
much smaller than the neutrino mass scales observed in the atmospheric and solar neutrino 
experiments. Furthermore, the upper bound on m ui becomes smaller for a larger M\. This 
clearly excludes the possibility that the three active neutrinos are nearly degenerate in mass, 
and indicates that one could find the absolute values of the heavier active neutrino masses 
from the results of neutrino oscillation experiments [1]. These conclusions remain valid even 
if there is entropy production as long as non-thermalized sterile neutrinos play the role of 
dark matter. 

8. Conclusions 

In this paper we have studied the production of sterile neutrinos in the z/MSM through active- 
sterile neutrino transitions from a thermal plasma, which does not contain any significant 
asymmetries related to lepton numbers. Though a first-principles formula for the production 
rate exists, which is valid to all orders in the strong coupling constant, its practical evaluation 
is subject to a number of uncertainties related to the strong interactions. The purpose of this 




(7.7) 
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paper has been to present a realistic "mean" evaluation of the hadronic contributions, and 
to estimate conservatively all the uncertainties that this mean result is subject to. 

The most important result of this paper is encoded in the four lines shown in Fig. [H They 
correspond to the case when there is no entropy production (S = 1) due to the decay of the 
heavier sterile neutrinos of the zaMSM. The area above the dotted line is certainly excluded: 
the amount of dark matter produced would lead to the overclosure of the Universe. The region 
below the dashed line is certainly allowed: the amount of sterile neutrinos produced due to 
active-sterile transitions is smaller than the amount of dark matter observed. Any point in 
the region between the two solid lines (corresponding to the "most reasonable" model for 
the hadronic contributions that we have been able to come up with) can lead to dark matter 
generation entirely due to active-sterile transitions. A maximal variation of the parameters 
of our hadronic model extends this region to the space between the dotted and dashed lines. 
In the case of entropy production with a factor S > 1.0, all these four lines simply move up 
by a factor S. 

As Fig. [9] shows, active-sterile transitions can account for all of dark matter only if M\ < 3.5 
keV, if the "most reasonable" hadronic model is taken. The most conservative upper limit 
would correspond to Mi < 6 keV, if all uncertainties are pushed in the same direction 
and also if the most stringent X-ray bounds are relaxed by a factor of two. Therefore, if 
the Lyman-a constraints from refs. \32\ [33] are taken for granted, the production of sterile 
neutrinos due to active-sterile neutrino transitions happens to be too small to account for the 
observed abundance of dark matter. In other words, physics beyond the i/MSM is likely to 
be required to produce dark matter sterile neutrinos. Another option is to assume that the 
Universe contained relatively large lepton asymmetries [8]. We would like to stress, though, 
that (apart from the astrophysical uncertainties related to the Lyman-a data) the simulations 
mentioned have not utilised the correct non-equilibrium momentum distribution functions as 
given in Fig. [7J and may thus contain systematic uncertainties. 
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Appendix A. Two-particle phase space integrals 

We discuss in this appendix how the phase space integrals of Eq. (|3.ip can be reduced to a 
one-dimensional integration in order to allow for a numerical evaluation. 

Let P be a time-like four-vector, P 2 > 0. We can then define a Lorentz-boost, Ap, which 
takes us to the rest-frame with respect to P: [ApP] 1 = 0, i = 1,2,3. We denote Q = ApQ, 
with the components Q = (q ,q). To be explicit, this Lorentz-boost is given by 



Q 



o 



7 [q° + f3e ■ q] , (A.l) 

q = q - e • qe + j[f3q° + e • q]e , (A. 2) 

where 

= - l 4, 7 = ^=L=, e=£- (A.3) 

p° - P 2 IpI 

Denoting by f[u; Pi; P% Q] a generic Lorentz-scalar, for instance 

f[u;P 1 ;P 2 ;Q]=2Q-P 1 n B (u-P 2 )[l-n F (u-P 1 )} , (A.4) 
as would be relevant for the second term in Eq. (|3.ip . we then consider the integral 

1 s Jw?k Jw?k (2 ' )4{<4,(P! " Pl " Q) /|!,; Pl • Q] • (A ' 5> 

where Q = (<7°,q), P% = (Ei,pi) are on-shell four-vectors. Making use of the Lorentz-boost 
introduced above, we rewrite the integration in a frame where Q is at rest, q = 0: 



/ d3 P\ / ^P 2 . (2vr) 4 ^ 3 )(p 2 - Pl )(5( J E 2 - E 1 - Mrf f[u;P 1 ;P 2 ;Q] 
J (2^)32^ J (2^)32^2 K J K 1J V >J[ ' ' J 



V'=A Q V 



(A.6) 



The integral over p 2 is now trivial: 

i r d 3 P i 



/ r% TT% Siy/pl+ml-y/pl + ml-M!) x 

' AJpf+mfJpf + m^ 



x f[A Q u; (ypf + mf, Ipil^pJ; (ypf + m^, |pi|0^); A Q Q]| . (A.7) 



dxx2g(x) ^(y^T^-y^T^-M) = g( m2 - mi -M) Pl2(M) ^ 12(M)) , 

x 2 + w x 2 + m| " 



Moreover the integral over |pi| can be performed by making use of 
,,,s 2 < 

/0 \J X 2 + 77^ ^ _, 

(A.8) 

where the function 

Pij(M) = -i-yV 4 - 2M 2 (mf + m]) + (m? - m 2 ) 2 (A.9) 



2S 



is real and positive for < M < \m\ — 777.2 1 and M > nii + nij. The integral thus becomes 

' PJ ' 6»(77l 2 - 7771 - Afj) X (A. 10) 



(4vr) 2 M 7 

x y"dfi pi f[A Q u; (ypfTmf, Ipil^pJ; (^/p? + m|, |pi|O pi ); (Mj, 0)] 



|pl|=pl2(M/) 



Now, the integrand in Eq. (jA.lOj) depends on two spatial vectors, q, f2 Pl , and thus on one 
angle. We can choose q = (0,0, |q|), S7 pi = (sin 6, 0, cos 9), JdO pi = 2ir J^d9 s'rnO. However, 
the Fermi-distributions in Eq. (|A.4h do depend on 6, because the spatial part of Aqu is non- 
zero and proportional to q, so that the remaining integration is non-trivial. On the contrary, 
scalar products such as Q ■ Pi, Q ■ P 2 , Pi • P2 are independent of 9, because we can evaluate 
them in the "hatted" frame where q = 0. 



Appendix B. Three-particle phase space integrals 

We discuss in this appendix how the phase space integrals of Eq. f|3.9j) can be reduced to a 
three-dimensional integration in order to allow for a numerical evaluation. 
Denoting by f[u; Pi; P 2 ; P3; Q] a generic Lorentz-scalar, for instance 

f[u;P 1 ;P 2 ;P 3 ;Q}=T l n F (u-P 1 )n F (u-P 2 )[l-n F (u-P 3 )] , (B.l) 

where % = T\ or T 2 from Eqs. (|3.7j) or (|3.8I) . we consider the integral 



1 s Iiwk Iiwk Iiwk (2 ' )45(4,(Pi + P2 - Ps - Q) /|a; '••='*«•• 

(B.2) 

where Q = (<7°,q), Pi = (Ei,pi) are on-shell four-vectors. Making use of the Lorentz-boost 
introduced in Appendix A, we rewrite a certain Lorentz-invariant subpart of the integration 
in a frame where P3 + Q is at rest, p3 + q = 0: 

d 3 P3 f f d 3 pi f d 3 p 2 



(2tt) 3 2P 3 U (2vr) 3 2Pi J (2vr) 3 2P 2 

x (2TT) 4 5^\p 1 + P2 )6(E 1 + E 2 -E 3 -q )f[u ] P 1] P 2] P 3] Q}\ . (B.3) 

J 1/eAp, . oV 



The integral over p2 is now trivial: 



x /[Ap 3 +Q77; (Jpl + m\, |pi|^ Pl );(VPi + ra 2 HPil^pi); A P3+Q- p 3;Ap 3+ Q(3] 



(B.4) 
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Moreover the integral over |px| can be performed by making use of 

dxx 2 g(x) f— 2 I 2 pi 2 ( J B)5'(/5i2(£ ; )) 



r V dxx g (x) (5(v / x 2 + m 2 + v / x 2 + m 2_^ ) = g(^_ mi _ m2 )/ 

o • 'x 2 + to 2 a/x 2 + m\ 



E 



(B.5) 



where the function is defined in Eq. (|A.9p . The integral thus becomes 

f d 3 p 3 f 1 /■ ,M) ^12(^3 + g°) v 

x / [Ap 3 +o u ; ( \/pf + m i 'I pi I n Pi ) ; ( \/pi + m 2 » - Ipi I ^pi ) ; a p 3 +q p 3 ; ^p 3 +qQ] } > (b.6) 

where |pi| = pi2(-E3 + 9°)- 

Now, the integrand in Eq. (jB.6[) depends on three spatial vectors, q, P3, f2 pi , and thus 
in general on three angles. However, the Fermi-distributions only depend on two angles, so 
that the dependence on the third angle is very simple and can be handled analytically. To 
implement this in practice we may for instance note that, due to 0(3) invariance, we can 
choose q= (0,0, |q|) and P3 in the (x, z)-plane, P3 = |p3|(sin#,0,cos#), so that 



roc r7T 

2ir / |p 3 | 2 d|p 3 | / sm6d8 . (B.7) 
Jo Jo 



Defining the unit vector e = (P3 + q)/|p3 + q| as before [it now lies within the (x,z)- 
plane], we parametrize the remaining vector £lp 1 by using spherical coordinates 6,4> with 
e as the polar axis. In the original frame, f2 pi is then given by Op 1 = (e z sin 6 cos <f> + 
e x cos 9, sin 9 sin <f), e z cos 9 — e x sin 9 cos 4>) . With this parametrization and given the form of 
/ that appears in Eq. (|B.1|) . it is not difficult to realise that the integrand depends on <f> only 
as a 2nd order polynomial in cos (ft: f{4>) = a + bcoscj) + ccos 2 0. Therefore we can replace 



2tt 







Only a three-dimensional integration (over | p>3 1 , 9, 9) needs hence to be carried out. 

Finally we remark that in the other channels (3 — ► 1, 1 — ► 3), the role of P3 + Q is played a 
difference, for instance P3 — Q. For arbitrary P3, Q this is not necessarily time- like. However 
the (5-functions appearing in these channels always restrict the differences to be equal to a 
sums, for instance P1+P2, which are time- like. Therefore non-zero contributions only emerge 
from regions of the phase space where P3 — Q is time-like, and the procedure described above 
can be taken over with minimal modifications. 
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